% TSP results
clear
load tsp_data


num_coef = 3;

tempindex = find(year >1970 & year < 1987);
cons = ones(length(tempindex),1);
[temp,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
arvals(1) = temp(2);
temp = sqrt(diag(inv(length(cons))*(r'*r)*inv([cons L_tsp(tempindex) year(tempindex) ]'*[cons L_tsp(tempindex) year(tempindex) ])));
arsevals(1) = temp(2);


tempindex = find(year < 1971);
cons = ones(length(tempindex),1);
[temp,~,r,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
arvals(2) = temp(2);
temp = sqrt(diag(inv(length(cons))*(r'*r)*inv([cons L_tsp(tempindex) year(tempindex) ]'*[cons L_tsp(tempindex) year(tempindex) ])));
arsevals(2) = temp(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = 7; beta = .95; M= 1;

R = length(arvals);
pxest = zeros(2,1);
gammatsp  = zeros(R,1);
dgammatspdrho = zeros(R,1);

for i = 1:R
    pxest(2) = arvals(i);
    [gammatsp(i), dgammatspdrho(i)] = gammafunc_se(T,beta,M,pxest);
end



gammatspse = arsevals'.*dgammatspdrho; 

alphaf = 18.23./gammatsp

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Robustness Check 1
tempindex = find(year >1970);
cons = ones(length(tempindex),1);
[temp,~,~,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) ]);
arval = temp(2);
pxest = zeros(2,1); pxest(2) = arval;
gammatsp_rob1 = gammafunc_se(T,beta,M,pxest);

% Robustness Check 2
load tsp_data_rob


tempindex = find(year >1970 & year < 1987);
cons = ones(length(tempindex),1);
[temp,~,~,~,~]= regress(tsp(tempindex), [cons L_tsp(tempindex) year(tempindex) LL_tsp(tempindex)]);
ar2val = temp;
pxest = zeros(4,1); 
pxest(2) = ar2val(2);
pxest(4) = ar2val(4);
gammatsp_rob2  = gammafuncAR2(T,beta,M,pxest);